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Abstract 

We propose in this paper an adaptive reduced order modelling technique based on domain 
partitioning for parametric problems of fracture. We show that coupling domain decomposition 
and projection-based model order reduction permits to focus the numerical effort where it is most 
needed: around the zones where damage propagates. No a priori knowledge of the damage pat- 
tern is required, the extraction of the corresponding spatial regions being based solely on algebra. 
The efficiency of the proposed approach is demonstrated numerically with an example relevant to 
<""j , engineering fracture. 

Keywords: model order reduction, proper orthogonal decomposition (POD), domain decompo- 
sition, nonlinear fracture mechanics, system approximation, parametric time-dependent problems 

1 Introduction 

Engineering problems are very often characterised by a large ratio between the scale of the structure 
and the scale at which the phenomena of interest need to be described. In fracture mechanics, the 
initiation and propagation of cracks is the result of localised microscopic phenomena. These phenomena 
are usually represented in an homogenised manner at a scale which is suitable for the simulation: 
the scale of the coarser material heterogeneities (meso-scale) , or the engineering scale when such a 
coarse representation allows for predictive results. In any case, the local nature of fracture leads to 
large numerical models because sharp local gradient need to be correctly represented or because the 
meso-structure needs to be described in an explicit manner. To some extend, the availability of super- 
computing facilities alleviate this difficulty. However, in engineering design processes, a prohibitively 
high number of solutions might be of interest, for a range of values of design parameters, or to take 
into account the effect of randomness in the model for instance. Therefore, one needs to devise efficient 
strategies for the solution to parametric multiscale problems. In doing so, the availability of a range 
of efficient numerical methods for the solution to one particular realisation of the parametric problem 
(homogenisation techniques, advanced discretisation tools, domain decomposition and multiscale-based 
preconditioners for parallel computing) should not be ignored. 

Model order reduction techniques that are based on the projection of fine scale problems in reduced 
spaces are a potential solution to this issue. Such strategies rely on the fact that the solutions to the fine- 
scale problem obtained for different values of the input parameters can be often represented accurately 
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in low-dimensional subspaces spanned by well-chosen basis functions at the fine scale. Applying this 
idea, the numerous unknowns that arise from the discretisation of the fine-scale problem are reduced 
to a few state variables (i.e.: the amplitude associated to each of the basis functions). Of course, 
obtaining the aforementioned global basis functions still requires heavy computations at the fine scale. 
Therefore, this class of methods is of interest if (i) the goal is to interact with a model (one can afford 
expensive "offline" computations in order to allow the user to interact with the reduced model in real 
or quasi-real time) or (ii) the cost of computing the global basis remains small when compared to the 
cost of solving the fine-scale problem for a large range of input parameters. This paper addresses the 
latter case, with a restriction to the design of structural components under extreme loading conditions. 

Projection-based reduction methods have been extensively studied in system engineering (see the 
review proposed in []J ) , fluid mechanics [H 0, 0, 0, 0] and structural dynamics 0, 0, 0[, E3, Hll.ll2l| . The 
theory and applicability of various projection-based model order reduction methods such as component 
mode synthesis [l3|, 0] , the reduced basis method [3, [l5[ , the proper orthogonal decomposition [lfl 
[l7l 0] which will be used in this work, the a priori hyperreduction method [3, [l9[ or the proper 
generalised decomposition 2(J 21 , 22 1 are now well-established in the linear to mildly nonlinear cases. 
Some attempts have been proposed to extend this concept to strong nonlinearities, in particular in 
structural mechanics 23, 3, 24, 2||. This background makes it conceivable to use such methods in 



complex engineering problems such as fracture mechanics. 

Fracture mechanics is characterised by an intrinsic lack of separation of scales between the engi- 
neering scale and the scale at which damage initiation is described. Consequently, these problems are 
not directly reducible by the aforementioned methods (this fact will be illustrated in the core of the 
paper). More precisely, the level of reducibility of such multiscale problems depends on the region of 
the domain which is considered. Typically, the solution in the zones where damage initiates and prop- 
agates will not be correctly approximated in low-dimensional subspaces. To circumvent this difficulty, 
the idea followed in this work is to use a partition of the structural components into substructures 
and perform a reduction of the resulting subproblems only if such a reduction can be done without 
sacrificing accuracy. 

The concept of local reduced basis itself is not new. It probably originates from the work of Craig 
and Bampton 0, who proposed a reduction by projection on a modal basis defined over predefined 
subdomains. This idea has been explored and improved in [2allll.ll 21 . or coupled with other reduction 
methods, as in the case of the proper generalised decomposition [2o| . A closely related family of solvers 
uses this concept within local/global approaches: only part of the domain is reduced (sufficiently far 
away from the sources of nonlinearity) [l(| H3, 28, f|, or the global reduced model is locally enriched 



by a fine-scale description |29|, [3CJ, |3l| (these two approaches are equivalent when the reduced model 



is used as a prcconditioncr for the local fine-scale problem in the former group of methods [28J]). The 
work presented here is novel in the sense that (i) it is the first formal coupling between Schur-based 
domain decomposition approaches and model order reduction by the Proper Orthogonal Decomposition 
and (ii) it is, to the authors' knowledge, the first application of systematic partitioned model order 
reduction for multiscale fracture. 

Reduced order models obtained by the proper orthogonal decomposition (see for instance 32, 3^, 



34 . |35| ) are powerful tools to reduce the computational burden associated with the repetitive 



analysis of parametrised nonlinear problems. The principle is to build the projection basis from the 
knowledge of a set of fine-scale solutions corresponding to a small number of chosen values of the 
input parameters (the so-called "snapshots"). The proper orthogonal decomposition (POD) is used 
to extract attractive reduced spaces from these fine scale solutions in an "offline" phase (we use here 
the terminology developed for interactivity). Classical projection-based reduction is finally deployed to 
compute solutions to the problem for arbitrary values of the input parameters at reduced cost ( "online" 
phase). Let us emphasize the fact that, by construction, this family of reduction techniques rely on 
the "offline" computation of fine-scale solutions (like the reduced-basis method, and as opposed to 
the proper generalised decomposition and a priori hyperreduction methods, which only require cheap 
fine-scale predictors). 

These "offline" computations are potentially expensive in the case of multiscale problems, and our 
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conception of the design process is that domain decomposition methods 3(1 37 , |39|, which are, 
to date, probably the most efficient family of parallel solvers, could be used to make them tractable. 
Examples of parallel computations using domain decomposition methods in the case of fracture can 
be found in (40l . l4l| . The purpose of this work is to reuse the substructured nature of the information 
generated during the "offline" stage to accelerate the solution process of the "online" stage. The choice 
of the domain decomposition method itself is not of prime interest here. Conceptually, we believe that 
the work presented in this paper can be extended to Schwartz-based methods, as done for the proper 
generalised decomposition in the LaTin framework [201 ] , or to other Schur-based domain decomposition 
methods such as FETI, as presented in 11| for component mode synthesis. We will focus in this work 



on the primal Schur-based domain decomposition method proposed in 3y, |37[. This method relies 
on a static condensation of the subproblems on the interface degrees of freedom, and a solution of 
the resulting problem by a projected, preconditioned conjugate gradient in order to ensure a certain 
level of scalability. We propose to use the snapshot POD method to construct reduced models of the 
sub-problems corresponding to the interior degrees of freedom of each subdomain. 

The proposed substructured approach to model order reduction (see a schematic representation in 
figure [IJ is adapted to the multiscale nature of fracture problems and provides benefits in terms of 
applicability of POD-based reduction techniques, along the following lines. Firstly, the POD transform, 
even when using the snapshot technique proposed in [2| can be prohibitively expensive to compute. This 
issue was treated in [3| by preserving the distributed nature of the snapshot data and reconstructing 
an approximation of the first modes of the global POD transform from local transforms computed 
independently for each subdomain. In our case, the POD bases will be used locally, and therefore, 
their parallel construction is natural. Secondly, using local reduced bases means that the dimension 
of the reduced spaces, can be adapted to the level of nonlincarity of the subproblems (seen as a 
statistic correlation of the snapshot data by the POD transform). As mentioned previously, the domain 
decomposition framework makes it natural to switch from a model order reduction solver to a full scale 
solver for the solution of subproblems for which no relevant low-dimensional reduced space can be 
constructed. Notice that similar ideas have been used in the context of domain decomposition methods 



without reduction for the treatment of localised nonlinearities arising in fracture mechanics. In 42j . 
subproblems corresponding to domains far away from the zones of interest are treated as linear, and 
the local finite element discretisation is coarsened to allow for computational savings. In 43J and 411 ] . 
the preconditioner of the domain decomposition method is used for the coarse solution to subproblems 
that are far away from the process zones. At last, we believe that the systematic decomposition of the 
domain makes the solution of propagating nonlinearities by reduced order techniques more amenable 
than local refinements around evolving zones of interest. 

The paper is organised as follows. In section [2j we give the general assumptions regarding the 
class of nonlinear problems which are addressed in this paper, section [3] introduces classical model 
order reduction by projection. We focus on the snapshot POD methodology and establish the state- 
of-the-art of system approximations for nonlinear problems. An example of application of POD-based 
model order reduction in the case of fracture mechanics is presented to highlight the difficulties due to 
the local lack of correlation in the data. In section 01 we introduce the primal domain decomposition 
method, and formally develop a POD-based model order reduction of the sub-problems in a Galerkin 
context. An inductive method is proposed to determine the set of fine-scale solutions that should be 
used to obtain a certain level of accuracy in the partitioned snapshot POD. A system approximation 
strategy for the partitioned POD approach is developed in [5] Finally, we propose results in terms of 
running time in section[5](as a first step, the partitioned POD is used in a serial computing approach), 
and discuss further improvements for the proposed strategy. 



2 General problem statement 

We consider the evolution of a structure described by the partial differential equations of continuum 
mechanics (mechanical equilibrium and constitutive law with appropriate boundary conditions) on a 
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Construction of partitioned reduced order model 



Compute particular realisations Partitioned reduced basis 




Solution for arbitrary parameter using reduced model 




no reduction 

Figure 1: Schematic representation of the partitioned POD-based model order reduction strategy. The 
Snapshot POD transform is performed locally for each subdomain. The resulting subspaces are used 
as a projection basis for model order reduction of the local problems in the balancing domain decom- 
position method. If the convergence of the local POD transforms is not satisfying, the corresponding 
subproblems are solved without reduction. 
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bounded spatial domain J7, over time interval T — [0,T]. The evolution in time is supposed to be 
quasi-static. We focus on nonlinear constitutive material models representing the progressive failure 
of structure, such as plasticity or damage. We assume that these processes are rate- independent. 
The mechanical problem is parametrised by a set of real variables fj, that evolves in parameter space 
V c R n " . ~~ 

Performing a space discretisation (finite clement in our examples) of such a problem leads to a 
system of coupled nonlinear (ordinary differential of viscous effects are described) equations. We look 
for the parametric evolution of the state variables U(t, /z) £ R™ u satisfying the following semi-discrete 
problem 

V(t,,i)eTxP Fint((U(r,At)) Te[0>t] ,/i)+F ext (i,M) = 0. (1) 

The vector of internal forces, F int £ IR™ U , is a non-linear function of the current state variables U(t, //) 
(e.g. vector of nodal values of the displacement field in finite element; we will therefore refer to it 
as "displacement"). n u is the number of spatial unknowns of ([T]). As we model structural damage, 
the vector of internal forces also depends on the history of the state variables (U(r, A*)) Te ro t{ over tne 
past time interval [0,t[. Typically, the dependence of F int to the history of the displacement is due 
to non-reversible material processes. In the context of parametric problems, F int may additionally 
depends on the design variables (design-dependant elastic constants for instance). F ext e K" u is the 
vector of external forces, which may depend on time and on the design variables (design-dependant 
external load for instance). 

A classical time discretisation of semi-discrete system ([TJ is performed. We search for a sequence 
of solutions (U_(t n , Li)) t e j-h, where we introduce the discrete time space T h = {to, t\, t„ t } such 
that to = and t nt = T, which satisfies the fully discrete set of equations 

V (t n , fi)eT h xV F int ((U(i m , /i)) meI0>n] , A') + E ext (i„, At) = (2) 

System ([2]) is solved sequentially in time, and we assume that the structure is undamaged and at 
rest to- At an arbitrary time step t n , the discrete history of the displacement (U(i m , M)) m e[o,n-i] i s 
known, which permits to calculate solution U(t n , li). For readability, the dependence of the system 
of equations to the discrete history of the variables, to the time step and to the parameter will be 
explicitly written only if necessary. 

Discrete system @ at current time step t n is nonlinear. It is solved by a usual Newton-Raphson 
algorithm. At iteration i + 1 of the nonlinear solver, a tangent linear system is solved: 

K* AU' +1 = -R* . (3) 



where AU =LP +i -lT is an increment in the displacement vector (U is the unknown solution) , 
is the tangent operator and R 1 = F.i nt (U l ) + F e xt ^ s the residual. The Newton 



dU 



u=u« 

algorithm is stopped if the relative euclidean norm of the residual at iteration i + fe — St, is lower 
than a chosen tolerance e new . In order to avoid heavy notations, we rewrite equation @ in the form 
KV = F. 



3 Model Order Reduction and Proper Orthogonal Decompo- 
sition 

Let us recall that our goal is to solve problem @ for a range of values of design parameters. In this 
context, the property underlying the applicability of POD-based MOR is that variations in the design 
variables generate variations in the solution which can be represented in an attractive low-dimensional 
subspace of R" u . Supposing that we can obtain a basis for this subspace, for instance by a particular 
application of the Proper Orthogonal Decomposition, then the evolution problem @ , for any value of 
the parameter, can be solved at reduced costs by looking for the solution in the reduced space. 
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3.1 Projection-based model order reduction 

Let us write that the solution to J2j can be approximated, at any time and for any value of the 
parameter, in a known subspace of M 71 " of small dimension spanned by a few basis vectors (CJjgnnj 

V(i„,/i)er k x? U(tn^)^^C 4 (t n , t j,)a i (t n ,fi)=Q(t n ,p,)a(t n ,^). (4) 

i=l 

where C(t n ,/i) £ K™ u x R" c is a matrix whose columns are the basis vectors (C i (i n ,/x))j 6 [i j „ c ] and 
a(t n , m) is a vector of reduced state variables (ai(t n , M))ie[i,« c ]- We emphasize that the reduced space 
Im(C) is not necessarily independent of time and parameter but, at this stage, is known. 

Injecting this approximation into © at a particular point (t n , fj,) of the time-parameter space, one 
obtains an over-constrained problem in the n c reduced state variables a (n c << n u ) 

Eint (Cg) +£ext — R(<*) ( n u equations, n c unknowns). (5) 

Problem ([3]) can be solved in different ways, depending on the physical quantities of interest and 
on computational tractability issues. The most used methods are the Galerkin framework and the 
least-square approach. The latter reads: 

a = argmin ||R(a*)|| e , (6) 

a*GM"<= = 



where ||R||e = yR T 8R denotes a 0-norm of the residual vector R (0 is a symmetric, positive 
definite operator). Alternatively, in a Galerkin projection framework, a is the solution of 

£ T R(a) = 0. (7) 

We use the Galerkin approach. Nonlinear problem ([7]) can be solved by a classical Newton algorithm. 



The linearisation of reduced problem ([7]) at iteration i + 1 of a Newton solver (see 24j for more details) 
leads to the following problem: 

(Find 5a \ C T (F - KCfo) = o) <^ 5a = argmin ||F - KCSa* II x , (8) 

where 5a = a l+1 — a 1 . Linearised system (|8|) is a Galerkin reduction (or a least-square reduction as 
these two approaches are equivalent for the linearised problem when using a K _1 -norm) of linearised 
equation © with the kinematic constraint V = C 5a. The solution to (jHJ reads 

= (Q T Kg) 1 £ 7 f - (9) 

providing the reduced linearised operator (of very small size n c ) is invertible. 

At this point, we can notice the two following classical issues in projection-based model order 
reduction: 

• The well-posedness of tangent problems © and the accuracy of the solution strongly depends 
on the choice of the reduced space. 

• The Galerkin projection framework presented previously is inefficient. The tangent and residual 
of the initial problem of evolution must be evaluated at each iteration of the Newton solver. The 
evaluation of nonlinear function F int requires a global integration over domain f2. As a result, 
the numerical complexity of the reduction technique does not depend on the dimension of the 
reduced space but on the size of the initial problem, which results in insignificant speed-up. 

Therefore, a reduction method should provide a "good" reduced space (in the sense of accuracy 
and stability of the solution), and an associated reduced basis, as well as an "efficient" strategy for 
the solution of (|5|) (significant speed-up compared to the full model, without sacrificing the accuracy 
expected when using a good reduced space). These two points are discussed in the following sections. 
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3.2 Proper Orthogonal Decomposition in projection-based model order re- 
duction 



3.2.1 Proper Orthogonal Decomposition 

The proper orthogonal decomposition (POD) is a popular transform which is classically used to gen- 
erate relevant bases for projection-based reduced order models. Applied to our parametric evolu- 
tion problem, the POD decomposes the solution of the problem over the full time-parameter space 
V = T h x V as 

V(t„,/i)6P, U(t„,/i) = U(t„, ) u) + e(t„,/x) (10) 
"i> 

H(i„,/i) = ^2^.ji(t n ,n) =£7(t„,/x) , 
j=i 

such that U is the function of separable form (jTU)) that is the closest to the exact solution, 

U= argmin d(U,U*) (11) 

S* e{z | z(t„,/i)=^T(tn,M),V(t„,Ai)e^} 

with the metric d defined on the space U of functions defined over V with values in K nu : 

d: UxU -> M_ 

(U,U) h> d(U,U) ( > 



rf(U,U) = 




((/>.)i 6 ji inp ] are "space" vectors that belong to R" u , while (7i)ie[i,n p ] are scalar functions of time and 
parameter. We emphasise here the fact that the spatial basis (p is not known a priori but is assumed to 

be independent on time and parameter (i.e.: we perform a separation of variables). The POD essentially 

delivers a decomposition of the exact solution U into bi-orthonormal modes ((</>.), 7, (f n ,/z) ) 

V — 1 /ie[l,n p ] 

of decreasing importance. The truncation of those modes at order n p provides the best representation 
of the solution with a basis of n p modes in the sense that the sum over time-parameter space of all 
distances between the exact solution and its n p -order approximation is minimised. Distance rf(U, U) 
is expected to decrease quickly with the order n p of the decomposition. 

We can immediately see the benefits of such a transform for the construction of reduced mod- 
els. It provides a reduced spatial space span ((0)ig[i,„ p ]) , in which the solution over V is optimally 
approximated at any order n p (in the sense of distance (|T3")0 . 



3.2.2 Snapshot POD 

The POD transform (|10I13|) requires the knowledge of the exact solution over P, which is not compat- 
ible with our desired usage. However, one can derive a similar transform that computes an optimal 
decomposition of the solution U over a discrete subset V s = T h x V s of V . 

V(i„,^)G^ s , U(t n ,fi) =V s {t n ,fi) + e s (t n ,n) (14) 

rip 

U S (£„,^) =5^^ f 7»(*n,M) =^7(*njM)> 
<=1 

such that TJ S is solution to the optimisation problem: 

U s = argmin d s (U,U*) (15) 

U*£{Z I Z(t„,/i)=<#)-y(t„^),V(t„,Ai)G-P s } 
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with d s the metric defines on the space U s of functions defined over V s with values in 

<F : U s xU s -> E 



with 



(U,U) h+ d(U,U) (16) 



d s (U,£ S )=E ]T ||U(i n ,/i)-U s (t n)/ti )||l (17) 



P s = {fii, ... ,/in ff } is a discrete subset of the parameter space V. (U(i n , ^))t n eT h nev s are particular 
solutions of problem ([2]) for some parameters n <E V s , called snapshot [2j, for which the solution to the 
discrete evolution problem is available. The snapshot POD metric (|17j) can be viewed as a quadrature 
rule for ([13]) . 

Optimal reduced spatial space span f (^^iep.npjj ancl corresponding scalar weighting functions 
(7<)ie[i,np] are S iven b y : 

• 0. is the eigenvector of the correlation operator H associated to its i th largest eigenvalue A^. H 
is defined by 

M = E E U(t„,/i)U(* n ,/i) T . (18) 
t^ev t„er h 

• V (tntfJi) G V s , 7i (tn , H) = £Wn , M) 

The truncation error of a POD transform of order n p is given by 

rf s (U,U s )= g A*, (19) 

i=n p +l 

where n s is the number of snapshot solutions, and therefore the maximum possible rank for operator 

H 

The eigenvalue decomposition of H is obtained at relatively cheap costs when nt x < n u by 
exploiting the discrete nature of the available information (which is essentially the idea proposed 
by Sirovich in |2j). One computes the singular value decomposition (SVD) of the snapshot opera- 
tor H = (U(ti,//i) V{t 2 ,Hi) ... U(i„ t ,/x n J). The SVD reads § = QJi^ T with Q and ^ uni- 
tary matrices and S a rectangular matrix with diagonal upper block. We then have H = S S T = 
QSW T WS T Q T = QSS r Q T , which is the eigenvalue decomposition of H and the eigenvalues 

are the squares of the singular values of S. The values of the weighting functions (7i)ie[i.n P ] over V s 
can be readily extracted from matrix W if necessary, but this information is not of particular interest 
the present context. 

3.2.3 Reduced space in POD-based model order reduction 

The snapshot POD provides an optimal decomposition of the solution in the discrete space V s . It can 
be truncated at an order n p < n s for which the normalised truncation error 



E * 



^nap = rf S (U,U S ) = ^ , (20) 



i=l 

is sufficiently low. 
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POD-based reduced order modelling proposes to simply discard functions (7t)ie[i,n ] (which are 
only defined for a discrete set of parameter values anyway) , and look for the solution of the evolution 
problem for any value of parameter /i £ P, in reduced space span((</> i )i 6 ji j „ p j), the amplitude associ- 
ated with the basis vectors being computed optimally by the projection technique given in section [3~T1 
POD-based reduced order modelling is therefore an a posteriori interpolation method for the snapshot 
POD. In this context, we will write in the following that the space vectors of the snapshot POD are 
used to define a reduced space for reduced order modelling (which is therefore independent on time 
and parameter): 

V(t n ^)er ll xP,V«e[l,n c ] C i (*n,At) = i K = n p ) (21) 



Remark: A solution over the initial time-parameter space V could be reconstructed by an explicit inter- 
polation of the functions (7i)ie i j n ] (i.e.: interpolation by an arbitrary polynomial basis or by Kriging 
for instance, as proposed in \3A . \AA]). which would lead to a decomposition of type (jlOp . However, 
such an explicit interpolation approach in V is suboptimal as the behaviour of the balance equations 
between the pre-computed snapshot solutions is unknown. In addition, and this is probably the most 
important argument, the Galerkin projection framework defined in section \3J\ permits to reuse the error 
estimates available in finite element schemes for the certification of the implicitly interpolated solution 
(see fl2 . \sdl . Qj 0/ for instance), at least in the linear case. The extension of this idea to nonlinear 
problems is currently an active area of research and will not be addressed in this contribution. 

An important point to emphasise is the requirement to perform cost-intensive simulations to com- 
pute the snapshots. We assume in this work that the initial problem of evolution involves a large 
number of degrees of freedom in space and time and requires high-performance computing for an exact 
solution to be at reach. In particular, this solution can be obtained efficiently on parallel architecture 
by using domain decomposition methods, which are, to date, probably the best parallel solvers for 
structural mechanics. This requirement will actually serve our needs in the case of fracture, as shown 
later on in this paper. 

3.3 System approximation 

As stated in section 13. 1[ an approximation of the fully discrete system of equations (JSJ) must be 
associated with the choice of the reduced space. In order to limit the computational expense due to 
the evaluation of the nonlinear functions F int , two families of strategics have been intensively studied 
in the literature. 



3.3.1 Linearisation 

The first family proposes to linearise 46, 9|, or perform a higher-order Taylor expansion 47, 23, 48| of 
the nonlinear terms in the system of equations. The reduced linearised operators can be computed once 
and for all "offline" and reused "online" in the Newton solver. Obviously, the validity of the Taylor 
expansions is local along the trajectory of the reduced state variables. The authors of 47] proposed 
an elegant "offline" linearisation of the nonlinear terms of the discrete set of equations that depends 
on the value of the reduced state variables. In the "online" phase, the nonlinear terms of the discrete 
set of equations are approximated as a weighted combination of the "offline" trajectory-dependent 
linearisations. 



3.3.2 Evaluation of nonlinear terms on reduced spatial domains 

The second family of system approximations proposes to only evaluate the nonlinear function at partic- 
ular points of the domain. In a first subset of these strategies, the nonlinear function is reconstructed 
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by interpolation on another POD basis ("gappy" technique) [H, 0, IHj]. The expansion of the 
nonlinear term reads: 

Pint (Ca(t n , /*)) = ft(t n , /i) = H£(t n , /*) • (22) 

i=l 

The columns of D G l" u x R nd are spatial functions corresponding to a truncated snapshot POD 
expansion of F int (Ca). Static coefficients (3 are found by expressing that at any point (t n ,fj,) £ V, 
the interpolation is optimal for a limited number n sa of spatial degrees of freedom: 

jSfo.AO = argrrrm||D/3*(t„, M ) -F int (Ca(i„ l(l )) || E (23) 

- (5* — - 

P is a boolean diagonal operator with very few non-zero values n sa > corresponding to the evaluation 



degrees of freedom of the spatial interpolation of the nonlinear term. ||X||p = y X T P X is the semi- 
norm associated with P. Substituting this approximation into the full system of equation ([5]), the 
following reduced problem is obtained: 

= (= T =i=) 1 U T E: Eint (G a) + Eext = Efe) , (24) 

where D T PD is assumed to be invertible. This system can then be solved in a least-square sense 
or by a Galerkin projection. Only a restriction to the evaluation degrees of freedom of the nonlinear 
function is calculated to evaluate the residual of the system, allowing to speed-up the solution process. 



The second subset of these strategies, collocation methods [50, |51|, |52j, do not reconstruct the 
nonlinear function over the domain. They propose instead to look for a solution that is optimal with 
respect to a few of the equations of the initial system ([5])- This can be expressed in a least-square 
approach: 

a — argmin||R(a*)||p , (25) 
< t ■ 

or in the Galerkin framework 

£ T £R(a) = 0. (26) 

The strategies proposed in the literature for this two subset of techniques differ in the way of 
building P, which requires a critical trade-off between optimality, stability and tractability. 

3.3.3 Chosen strategy 

We will focus in this work on the "gappy" technique, as used in [35| and (4^|. Since the main objective 
of this paper is not the system approximation strategy but the introduction of the partitioned POD 
technique, this method is selected as the most widely used. We note for the following developments 
that a Newton iteration applied to the reduced system reads (in the Galerkin framework): 

Sa = (c 7 D(D T PD) D'PKC) 1 C T F ext , (27) 



This technique allows one to look for a solution in a reduced space at reduced integration costs. In 
, P is constructed such that the condition number of the reduced linearised operators is minimised. 
In the hyperreduction method [HH , the non-zero entries of P correspond to the largest values of the 
residual vector R(a) . In [49| , the points are selected to limit the growth of the residual error between 
a solution and its snapshot reconstruction. This is the technique we will use and it will be addressed 
in the last section of this paper. Meanwhile, we focus on the issue of computing and using relevant 
POD-based reduced spaces in the particular case of fracture mechanics, using a Partitioned POD 
approach. 
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3.4 Example of application of the POD in fracture mechanics 
3.4.1 Lattice model 

We consider a lattice structure made of damageable bars in uniaxial tension or compression. A bar 
marked b occupies a ID linear domain fiW embedded in R 2 , such that SI*- 6 -* = ft. n\, is the 

be[l,n b ] 

number of bars comprising the lattice structure. We will denote by {-Pi}ie[i,n pt ] the set of nodes of 

the lattice structures. Let us define the unit vector rS b ^ attached to bar b such that if Pi and Pj are 

p p 

the two extremities of QW and i < j, then = \\p.p.\\ ■ We denote the local coordinate of point 

M G by s<® = \\ PjM \\ . 

We look for a two dimensional displacement field q, and a scalar stress field N over fi that satisfies 
the system of equations given below. The restriction of these fields to bar b will be denoted by q^ 
and respectively. 

The local mechanical equilibrium equation of the lattice structure reads, at any point of fl^ b ': 

dN (b ) /i\ 

=0 (28) 

At a lattice node Pi between a set of bars S n ,ij the tractions must satisfy the equilibrium 

z^sg =o < (29) 

or 

£ A^lg+iVdip^O, (30) 

if Pi G where <9f2/ is the part of the boundary of Q where Neumann boundary conditions are 

applied. are tractions that are applied on dftf. In equilibrium equation (|29j) . nip, = —n^ if 

s[p. = (first extremity of the bar), and = otherwise (second extremity of the bar). At node 
Pi , the continuity of the displacement field must also be ensured 

£|A=£|fl V(M')eB M . (31) 

The displacement field satisfies Dirichlet boundary conditions on boundary d£l u = dft\dfif, which 
reads 

2|A=£d| fl VP* GSflt,. (32) 
The constitutive law, which is detailed below, reads at time t: 



where the deformation eW is defined by the compatibility equation 



3.4.2 Damage model 



The fracture of the lattice structure is described by classical damage mechanics [53j . We postulate the 
existence of a free Helmholtz energy 

iP(e,d) = ±E(l-d)Se 2 (35) 
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E is the Young's modulus of bar b, S is its section (assumed constant), and d is a damage variable 
that ranges from (safe material), to 1 (completely damaged material point). The state equations are 
obtained by derivation of the free energy with respect to the state variables. 

N= =E(l-d)Se, (36) 
ot 

Y is a driving force associated with the damage variable d. To close the system, a simple evolution 
law is formulated: 



d = minjmaxj^-j (38) 

Notice that the history dependency in the previous equation (non-reversibility of the damage process) 
is inherited by the discretised system of equations. 

The damage problem is solved by linear finite elements. The displacement field q is searched for as 
a linear combination of compactly supported hat "functions" . The nodes of the finite element mesh 
coincide with the nodes of the lattice structure. Substituting this approximation into the weak form of 
the balance equations, one obtains a semi-discrete system , which is then discretised in time to obtain 
©■ 



3.4.3 Parametrised fracture problem 

The leftmost part of the structure is fixed (null Dirichlet boundary conditions) while a prescribed 
displacement, which puts the structure in tension, is gradually applied on the rightmost part. The 
direction of the load is controlled by an input parameter 9 which ranges in V = [15°, 45°]. An initial 
crack (notch) is defined at the top centre of the structure by initially setting the damage variables of 
the corresponding bar elements to 1 , as illustrated in figure [2l As the load is progressively applied to 
the damageable structure, the crack propagates. The time evolution of the crack propagation problem 
is discretised using 10 homogeneous load steps. 



Initial crack 




6 € [15°, 45°] 



Figure 2: Definition of the nonlinear lattice problem used for the numerical experiments of this paper. 
The loss of stiffness of each bar while increasing local strain is described by a damage model. The 
direction of the prescribed displacement on the right-hand edge of the rectangular lattice structure is a 
parameter of the model. We want to be able to predict the propagation of the damage onset (initially 
damaged bars represented in black) for any angle of the prescribed load. 



Our goal is to predict the damage state in the lattice for any arbitrary angle 6 € [15°, 45°] without 
solving the full model for each angle. The solution for an arbitrary angle 9 will be sought after in a 
space generated by the spectral analysis (SVD) of precomputcd solutions for the two extrema values 
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of 9, 15° and 45°, obtained by solving the full, non-reduced, model. To be precise, 20 solution vectors, 
obtained for the first 10 time steps for each of the two loading angles, are used to build the snapshots. 
The normalised truncation error f snap of the snapshot POD as given in equation (1201) is set equal to 
10 -2 (see figure [3]), which leads to three orthonormal displacement basis vectors. Results are shown 
in figure HI 




Figure 3: Convergence of the normalised POD error indicator as a function of the order of truncation, 
for increasing size of the snapshot. The lack of correlation due to the crack propagation introduces a 
local error of projection. This error can be seen as a lack of correlation below a certain value of the 
snapshot POD error indicator. This threshold is relatively low due to the global nature of the metric 
used to evaluate the accuracy of the projection. 

It is noticed that each load angle 9 leads to a crack/damage zone propagating approximately 
orthogonally to the load direction, as is commonly observed in fracture mechanics. Consequently, 
each and every load angle leads to a different crack which cannot be well represented by a linear 
combination of the cracks obtained for a limited number of snapshot solutions (figure El bottom). 
In fact, the solution to parametric problems involving the evolution of topological changes cannot, 
in general, be obtained efficiently using a method based on the separation of variables (unless when 
manages to map the physical space to a reference space were correlation in the data can be retrieved) . 
However, the topological changes are localised in the case of fracture. In the regions that are far away 
from the crack, the solution is indeed well approximated by a linear combination of the pre-computed 
basis vectors. Consequently, a model reduction can still be performed but only over selected regions 
of the domain. 

The following section presents a possible strategy to implement this idea based on a domain decom- 
position method where the subdomains are selectively and independently reduced, based on a criterion 
described in section 1431 

Remark: The initial crack is meant to provide a stress concentration zone from which fracture will 
initiates. We emphasize here that this is an idealisation of a general situation in realistic engineering 
components. Cracks initiates from joints, supports, free edges (large shear stresses due to a mismatch 
between elastic properties in composite laminates for instance), non- smooth parts of the boundary of the 
component (corner), or from interior regions which are subjected to extreme stress concentration under 
particular external loading conditions. Therefore, the regions of potential initiations are not arbitrary 
for a given parametric problem. In the particular example treated in this paper, fracture propagates 
from the notch which was introduced in the geometry. However, in all the following developments, we 



13 



Construction of the reduced order model (ROM) 



Compute! particular 
realisations (snapshots) 

15" , 10 th (last) time step 



Reduced basis 



POO 



45° , 10"' time Step 



\ 
\ 
\ 
\ 




Solution at arbitrary angle using the reduced model 

1 



$ = 30" , m ili time step 



Solution of the ROM 



Solution to the 
full, unreduced, 
model 




Error 



Figure 4: Schematic representation of the Snapshot POD model order reduction technique for the 
proposed parametrised problem of fracture. The "Exact" time evolution of the problem is computed 
for several particular values of the parameter. A reduced space is extracted from this snapshot by a 
spectral analysis (SVD). In a second phase, the initial problem is at reduced cost by projection in this 
reduced space for arbitrary values of the parameter of the problem. In the case of fracture mechanics, 
the projection error localises in the "process zone" surrounding the crack. Far away from it, a reduced 
space of small dimension is sufficient to capture the solution with a high level of accuracy. 
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do not make use of the knowledge of the position of this initial defect, which emulates the existence of 
a priori unknown zones of stress concentration in the structure. 



4 Partitioned model order reduction approach 



4.1 Principle of the primal Schur-based domain decomposition method 

Schur-based non-overlapping domain decomposition methods (see a review in [54j]) are dedicated to 
the solution of large linear systems. In our case, we use the primal Schur-based domain decomposition 



(balancing domain decomposition (BDD) [36|, |37 



38[)) to calculate successive Newton iterates. Schur- 
based domain decomposition methods propose to condense the balance equations on the interface 
degrees of freedom (degrees of freedom that are shared by at least two subdomains), by eliminating 
the interior degrees of freedom. The resulting problem (called "interface problem" ) is of much smaller- 
size than the initial problem. It is solved by iterative solvers, usually preconditioned Krylov subspace 
algorithms, which are particularly well-suited to parallel computing. The derivation of the precon- 
ditioner for the iterative solution of the interface problem is a key point to obtain an efficient and 
scalable domain decomposition method. 




Figure 5: Subdivision of the domain of interest into 10 non-overlapping subdomains. Vj 6 ^ is the 
restriction of a vector V*- 6 -* of nodal values of subdomain e to the internal degrees of freedom of the 
subdomain, while Vv, corresponds to the interface nodal values. The superscript between brackets 
indicates the number of the subdomain. 
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Let us now give an overview of the BDD method. We look for the solution of a Newton iteration 
at time step t n and parameter [ij G V s , during the calculation of the Snapshots required to perform 
a POD transform. The discrete linear system to solve at the i th iteration of the Newton solver initial 
problem ([2]) is 

KY F (39) 
where the notations have been defined in section [2j Domain is now split into non-overlapping 



subdomains (Q^) 



ee[l,n c I 



such that (J 



e£[l ,n e 



£}( e ) = 57, as illustrated in Figure [5l We later refer 
to the subdomain set [1, n c J as £ for simplicity. Let \} e > be the vector of local degrees of freedom 
corresponding to fl.( e \ and N^ 6 - 1 the boolean operator that maps global degrees of freedom V into local 

degrees of freedom V (e) by V (e) = ^ (e)T V. 
The local equilibrium on reads 



(40) 



where K^ 6 -* = N*- 6 -* KN' e ' is obtained by a standard sub-assembly process (K^'s size is typically n c 
times smaller than K when the subdomains have all equal sizes), and vector A' 6 ' corresponds to the 
contact forces applied on the boundary dQ^ of by the adjacent subdomains. Only the entries of 
corresponding to dCl^ are non-zero. 
By separation of the interior local degrees of freedom Vj (of size n[ e ^ ) and interface local degrees 



of freedom V b (of size n^ 1 ) in the form 




the local systems can be written as: 



K (e) K (e) 
=11 ^=ib 



v (e) 



F 



(e) 



(41) 



(42) 



=bi =bb. 

The interior degrees of freedom Vj are eliminated from the local system by static condensation: 

v[ e ) = (f. (6) - Y b e) ) . (43) 

This leads to the condensed local problem 

F (e) , 



S (e) V (e) 

=11 D 



(44) 



where the primal Schur complement is given by = — K^K^ , and the condensed 

=p ^ =p =bb =bl =n =ib 

forces Fi e) are given by £( e) = F^ e) - K^K^^fH. 

— c o J — c — b = bi =ii — i 

In the BDD method, unique displacement interface degrees of freedom are introduced, such that the 
kinematic continuity is automatically ensured. It is related to the local displacement on the interface 
by the following equation: 

(45) 



A< e >V 



V 



(e) 
b ' 



where is a boolean assembly operator. 



Interface forces A b satisfy the interface equilibrium equations 

E4 (e) A b e) = o. 



(46) 
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The global Schur complement system is found by left multiplying each of the equations ()44j) by the 
corresponding assembly operator and summing up over all the subdomains: 

EA (e) ^ e) Yi: u E4 (e) ^ ) + E4 (e) ^ e) ( 4? ) 

ee£ ee£ ee£ 

= from l|46j 

^E (A {e % e) A {e)T ) Y b = E A (e) Fi e) • (48) 

e££ e££ 

In a compact form, we can write 



'S =yA( e »S< e »A l 
= P z_^ = =p = 



(e) T 



S V b = with { . (49) 

Global interface problem (j49|) can be solved in parallel iteratively using a Krylov-subspace method 
such as the conjugate gradient or GMRes (55|(or BiCGStab [H^|) in a non-symmetric case. In this 
framework, the global Schur complement need not be assembled. Instead, whenever it is needed in 
a matrix/vector multiplication, the multiplication is performed locally on each subdomain using the 
local Schur complements. The outcome of these local multiplications is then assembled: 

For any X b : SX b = E A (e) i|, e) A (e)T X b ■ (50) 

The local stiffness matrices K^f-', involved in the computation of the local Schur complements 
are inversed directly. Using this method it is possible to perform the matrix/vector multiplications 
(computationally the most demanding part of a Krylov-subspace method) in parallel, thus providing a 
computationally efficient strategy. In a similar way, the dot products involved in the iterative algorithm 
can also be performed in parallel. 

4.2 Formulation of reduced order modelling in the domain decomposition 
framework 

4.2.1 Local snapshot POD reduced spaces 

We propose to use POD-based model order reduction on the interior degrees of freedom of each 
subdomain. We assume that a snapshot {U(t, A*)}(t ^e-p 8 ^ s available. Local POD spatial bases 

( Ctt ) ,, of dimensions ni are computed for the interior degrees of freedom of each subdomain 

V '*/<e[i,n£°] 

e €E [1, n e ] as described in section [3] Accordingly, the normalised truncation error of the local snapshot 
POD transforms are defined as follows: 



E E 

\ snap J 



u i (/w')-£(a ( :fiJi(w))c l (e) 



E E llHi(in,|U)ll 



(51) 



where a partitioning into interior and interface degrees of freedom similar to (141 [) has been performed 
on the snapshot solutions {U_(t, /u)}(t M )e-p=- We define the local operator Cj e ) whose columns are the 
local POD basis vectors of subdomain e. 
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4.2.2 Local projection 

We look for the interior degrees of freedom in the local reduced spaces as follows: 



V 



(e) 



v (e) " 




"c.< e W e) " 


_ Y (e) 







(52) 



where Sa^ is the vector of local reduced state variables. Thus, we have to solve the following system 
of equations for each subdomain e e £: 



Kl e) K.t 3 ' 

^=11 ^=ib 

K (e) K (e) 

.^bi =bb. 



C (e) <5 



a 



p(e) 
T?(e) , \(e) 



(53) 



This is a set of n.; + equations in n[ e ' + unknowns and is thus overdetermined, as nf + 



,( e ) 



Consistently with the developments proposed in section [3~T| we perform a Galerkin 
projection: the residual of local system (|53p is required to be orthogonal to the reduced space, which 
reads 



C(e) 





T (S) 

=d 



p(e) , \(e) 

£-b +A b 



^=ii ^=ib 

K (e) K (e) 

=bi =bb 



V(e) 



0. 



(54) 



We end up with the following linear, symmetric system for the expression of the reduced local equilibria: 



o 

y(<0 



X (e) 



V 



(< : '> 



where < 



=i ^ii =i ^ib 

K (e) 
=bb 



Fi e ) 



n(e) 



(55) 



which is now a well-defined system of much smaller dimension (ni^ + n^) than the original one 
(n^ + n b ) in equation (|42)) . 



4.2.3 Condensed interface problem 

Similarly as described in section 14. 1[ local systems (|55|) are condensed on the interface degrees of 
freedom, and formally assembled. To do so, the reduced state variables Scx'f'' are eliminated: 



5a s (e) =K(f) -1 (F- (e) -K. ( f 3 V b e 

1 =n,r \ — =ib,r — D 



(56) 



where K.. = C^ T C[ e \ K., = C^ T K^. The resulting assembled, condensed reduced system 

^=n, r ^=11 ^=ib,r ^=ib 

reads: 



S V b =F C 

=p,r — u c 



with 



=p,r Z— < = =p,r = 

Fc, r =E= (e) F^ 



(57) 



with the expression of the local condensed operators = K.( e ) — K.^) K.^f' , the local con- 

— " — -ab ^=bi,r^=n,r ^=ib,r 



=p,r 



densed forces = F^ e) - K[ E) 'f.^ and K w . = K£?C. (e) . Problem (£7b can be solved using 

— ' — D ^bi,r^n,r — i)^ ^bi,r ^bi ^1 

a Krylov algorithm, as described in section I4TT1 

We can now go one step further and choose not to reduce the local problems corresponding to some 
of the subdomains. Indeed, if localised non-linearities arise (damage in our case), the local reduction 
based on the separation of variables might be inefficient: a prohibitively large number of spatial basis 
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vectors might be required to obtained the desired accuracy over the whole parameter space (recall the 
results of section l3.4.3[) . This particular issue will be addressed in section FOl So far, we will assume 
that the subdomains are divided into two complementary sets £ rcd U £ nrcd = £, where £ rcd is a set 
of subdomains for which reduction is numerically efficient, while £ nrod is the complementary set of 
subdomains, for which a direct solution to the corresponding local problem is preferred. The resulting 
hybrid condensed reduced problem formally reads: 

fs = J2 A (e) i (e) A (e)T + £ A (e) i (e) A (e)T 

l ee£ ,cd ee£=«= d 

4.3 Local error estimation by Cross- Validation 
4.3.1 Principle 

The partitioned projection approach described in section 14.21 allows us to construct reduced spaces 
that arc independent for each subdomain. We propose here a simple scheme in order to (i) determine 
independently the dimension of the local reduced space that is necessary to achieve a predefined 
accuracy for each of the subdomain (ii) evaluate whether a subproblem is reducible or not in the sense 
of the usual separation of variables assumed by the POD transform. 

These two points are addressed while considering that a relevant snapshot is a priori available. This 
relevant snapshot should explore the parameter space sufficiently. At the same time, one does not want 
to compute too many snapshot solutions, in order for the "offline/online" strategy to remain efficient. 
Ultimately, a third point has to be added for the design of a substructured learning strategy: (hi) assess 
whether the snapshot contains a sufficient quantity of information, and generate additional, well-chosen 
data if required. This last issue is extremely complicated to address. Some recent propositions have 



been made in [57|, |4J], but most of the studies on the POD, or the Principal Component Analysis in 
the statistics community (a recent review is provided in [Hij ]) consider that a sufficiently rich snapshot 
is available, and perform the spectral analysis without considering the need, or the possibility, to 
regenerate data a posteriori. 

We will here address points (i) and (ii), while point (iii) will be left to the perspectives of this 
work. The particular technique used in this paper relies heavily on cross-validation (CV, see [Hij ] in 
the context of the PC A), and more precisely the Leave-One-Out (LOOCV) technique. In order to 
validate the predictivity of statistical models, one usually divide the available data into a training 
set and a validation set. In our application, the training set is the snapshot: the set of solutions to 
the parametric problem of evolution that corresponds to parameter values in V s . The relevancy of 
the reduced spaces generated by the snapshot-POD can then be evaluated on a set of additional fine- 
scale solutions: the training set. Using independent training and validation sets permits to avoid the 
overfitting behaviour (or "Type-Ill' error" in statistics) that is classically observed in any regression- 
type model. In our context, the snapshot POD transforms only minimise the mean square error of 
projection of the snapshot solutions (fTTl) . Therefore, the associated error estimate (TlT)l) is expected to 
underestimate the error of projection associated to a hierarchically enriched snapshot, and in the limit, 
to underestimate the integral form (|13|) of the error of projection. Using different set of solutions to 
identify the reduced space and to compute the error of projection permits to avoid this effect, but at 
the cost of additional data, which means further cost-intensive fine-scale solutions in our case. 

The cross-validation error estimate avoids these additional computations by emulating the inde- 
pendence of training and validation sets using the same dataset. In order to do so, the summand in 
equation (f5Tj) is calculated using the local reduced basis obtained by a snapshot POD transform of 
all the available snapshot solutions but the one corresponding to the value of the summation vari- 
able. This is the usual LOOCV strategy applied to the POD. This can be written formally, for any 
subdomain e G £: 
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snap 



E E 



u,(t«,A*)-x;(s 



HiCw) Ci, 



E E iiu ; (^^)ii 



(59) 



vector (c[ '''I 



which is parametrised by the summation variable 



the modified reduced basis 

fi£P s , are the ric 6 " 1 first eigenvectors of the following modified POD operator: 

= ^ ^ U(t„,/x*)U(t n) /x*) T - 



(60) 



Technically speaking, the computation of this estimate requires to perform an SVD for each of the 
snapshot solutions (and for each subdomain). 

Let us remark that statistical error estimates are commonly used in the context of deterministic 
parametric problem. For instance, classical Kriging interpolations are based on a randomisation of 
the field to interpolate. We refer to 34j, |44j for recent combinations of Kriging and POD. The later 
contribution uses the LOOCV both as an error estimate and as a criterion to refine the snapshot space 
in a hierarchical manner. 



4.3.2 Application 

The LOOCV error estimate is now applied to the problem of fracture. The parameter space is sampled 
using a regular grid of 5 parameter values including the extremities of V = [15° 45°], which is, for 
the moment, assumed to be sufficiently fine for our purpose. In figure [6j the corresponding LOOCV 
estimate is plotted as a function of the dimension of the local reduced spaces for 4 different subdomains: 
subdomain 6, which is the most affected by the damage propagation, subdomain 4, which contains the 
"tip of the crack" for a range of parameter angles, and subdomains 2 and 7, which are further away from 
the source of nonlinearity (or lack of correlation, depending on the point of view) . Again, we emphasize 
that we treat all subdomains in the same manner. We do not make use of an a priori knowledge of 
the spatial distribution of damage. The lack of reducibility of certain parametric subproblems must 
be an output of the method. 

The effect of the localised damage on the error estimates of each subdomain is relatively clear. 
For subdomains that are far away from the crack, we observe a fast convergence of the LOOCV error 
estimate with the dimension of the local POD reduced spaces. A satisfyingly level of predictivity, set 
here to the threshold Psnap < 10 -3 , is obtained with 4 to 5 reduced basis vectors. It is interesting to 
notice that we do not obtain a clear "elbow" in the convergence curve, which is often used to define 
the "dimensionality" of the underlying parametric problem. This is, to our best knowledge, due to the 
far effect of the crack. The lack of correlation due to the local damage tends to pollute the remote 
area. Further evidence of this fact can be found in our recent investigations about this particular effect 
[60]. For the subdomains that contain most of the damage, the observed convergence curves are much 
slower. The required accuracy for subdomain 4 is obtained with 7 local POD basis vectors. In the case 
of subdomain 6, the LOOCV error estimate does not reach the predefined threshold. This indicates 
that the corresponding subproblem should not be reduced. 

We have now achieved our objective of choosing the dimension of the local reduced spaces based 
on a CV error estimate, and identifying non-reducible subproblems, based on an assumed sufficiently 
fine sampling of the parameter space. The local reduced spaces obtained in this section will be the one 
used in the following to demonstrate the numerical efficiency of the partitioned model order reduction 
approach. 
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Figure 6: Cross-validation error estimate as a function of the order of the POD transforms for 4 of 
the 10 subdomains. The snapshot comprises 5 instances of the solution to the parametric problem of 
evolution. 

5 System approximation in the partitioned model order re- 
duction approach 

5.1 Local "gappy" approximations 

We propose here to extend the concept of "system approximation" to the partitioned model order 
reduction introduced in section 2) As mentioned previously, we choose to apply a tailored version 



of the "gappy" reconstruction technique presented in different contexts in [5|, |49|, |35[, called (Dis- 
crete) Empirical Interpolation Method (DEIM). It is important to realise that the gappy technique 
approximates the Galcrkin projection framework described in section 2J Therefore, the system approx- 
imation will systematically be compared, or optimised, with respect to this reference framework and 
not with respect to the full-scale solution. This approach to system approximations is characterised 
as "consistent" in [35[ • The starting point of the gappy technique is to compute a local reduced basis 

|p_| e ' I e G £ rcd | to approximate the local vector of internal forces |Fj^ ; | e G £ rcd j, as detailed in 
the non-partitioned case in section [3j 

V e G £ Icd , Va G R™^ , p£> (g^ e) a\ « n[ e) ^ (a) (61) 



The static coefficients are obtained by minimisation of a distance between the previous approximation 
and the exact local vector of internal forces, evaluated at an arbitrary point along the reduced kinematic 
trajectory. This distance is measured at a set of sample spatial points only, which yields the gappy 
approximation: 

The local boolean operator Pj^ operating on the subdomain e G £ rcd is such that only the diagonal 
entries that correspond to all the degrees of freedom of a small set of internal nodes of e are set to 
one. These nodes are called "control points" or "control nodes" . We define the local gappy operator 

of subdomain e by G^ e) = pj e) (D^pWD^-iD^pW. 



21 



Upon linearisation of the local nonlinear subproblems (i.e.: derivation of the vector of internal 
forces with respect to the reduced state variables and interface degrees of freedom), and taking into 
account the gappy approximation (|62[) . one gets a modified expression of the local tangent systems 
(compare equation (|53p ): 



Q(e) j£(e) q(e) j^-(e) 
=DI =<bb 



(e) 
b,sa 



v 



(e) F (e) / r (e) (e),i 



-i,sa 

(e) 



) 



r (e) ' 
-ext,i 



b,sa 



(63) 



Recall that consistently with the notations introduced in sections [5] and SI F^ is minus the contribution 



of subdomain e to the restriction of residual F-^fU^' 1 ) + F^ t to the interface degrees of freedom. 
Exponent i denotes the previous iteration of the Newton algorithm. Index sa denotes variables that are 
obtained approximately using the system approximation, as opposed to the ones obtained in section U 
without the gappy technique. 

As mentioned in section [31 this overdetermined system can be solved in different ways. We use a 
Galerkin orthogonality condition, which, together with the gappy approximation, yields the following 
Petrov-Galerkin formulation for the tangent subproblem corresponding to subdomain e G £ rcd : 



.(<0 

-b,sa 



X (e) ' 

da) ' 

i,sa 

V (e) 
— b,sa 



K 



with < 



^ii =i ^ib 

j^(e)Q(e) j^(e) 

=bi =i =bb 



(64) 



A condensed interface problem can finally be obtained as follows: 



S V b =E C 
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The method to obtain the expression of modified primal Schur complement S ga and the correspond- 
ing condensed right-hand side is not detailed for the sake of concision. It follows exactly the method 
deployed to get their counterparts when no system approximation was used (see equation (|57I0 . 

Notice that the symmetry of the condensed interface problem is lost when using the gappy tech- 
nique. This issue can be alleviated by using a GMRes algorithm. 

The key benefit in using the gappy technique is that only the components of the local tangents 

and local residuals that are not filtered out by operators {pj e) |ee £ rcd } need to be computed. More 

precisely, the assembly of the tangents and residuals is performed via loops over all elements. With 
the system approximation, only contributions from elements that are connected to one of the "control 
nodes" are computed, which results in significant computational savings, as will be shown numerically 
later on. The set of elements over which an integration of the internal forces is required is called the 
reduced integration domain. An example of such a domain is shown in Figure [7] 



5.2 Construction of the system approximation 
5.2.1 Static POD bases 

To generate the local bases |^D^ e - ) | e S £ rod j , we develop a technique that is strongly inspired by the one 

proposed in [35j j . Equation (|61[) indicates that we would like the system approximation to be accurate 
for any set of local reduced state variables. However, we can reasonable restrict ourselves to the state 
variables that are observed on a set of particular solutions to the reduced parametric problem (without 
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Figure 7: Example of a reduced integration domain. Subdomain 6 is not reduced. Therefore, all the 
associated elements belong to the integration domain. Since the interface between substructures is not 
reduced in the proposed primal version of the Schur-based partitioned model order reduction method, 
all the elements that are connected to the interface also belong to the reduced integration domain. 
The remaining controlled nodes are obtained by a partitioned discrete empirical interpolation method. 



system approximation). In order to do so, one first solves the time evolution problem corresponding 
to snapshot space V s using the Galerkin framework described in section @] We emphasise the fact 
that the local solutions that are obtained in this fashion belong to the local POD reduced spaces. As 
mentioned in the previous subsection, this is the reference for the system approximation. We now 
want to approximate the spaces spanned by the local vectors of internal forces corresponding to the 
successive iterations of the Newton algorithm used to compute these reduced solutions. We call this 
space the "static snapshot". It can be represented mathematically, for any subdomain e € £ rod , by 
the following set: 



(66) 



In the previous expression, rinew denotes the number of iterations of the Newton algorithm used 
to solve the problem of evolution for parameter fi at time step t n . A singular value decomposition 
can now be used to compress and hierarchically order the information contained in this set, which 
is similar to the technique used to obtain the reduced bases for the displacements and constitutes a 
keystone for the greedy selection of the reduced integration domain proposed in [49[ . Technically, 
for each subdomain e G £ rcd , an operator whose columns are the vectors of set (po) is created, and is 
decomposed by singular value decomposition. The left-singular vectors associated to singular values 
that are larger than a certain tolerance define the columns of operator pj 6 * 1 (see below) . 



5.2.2 Dimension of the local POD "static spaces" 

One question that now arises is how to choose the order of truncation of the local SVD performed 
in the set of vectors of internal forces. In other words, one needs to choose the rank of the matrix 
of left singular vectors Dj e ' for each subdomain e G £ rcd . One could simply truncate the local POD 
expansions such that the truncation error becomes smaller than a predefined tolerance, as proposed in 
section T4.3I when defining the dimension of the local reduced spaces for the displacements. However, 
we prefer here to link the error generated by the gappy approximation to an error measured in terms 
of displacements, such that it can be compared to the error introduced by the truncation of the local 
snapshot POD performed to generate the local reduced spaces. 

In order to detail the numerical implementation of this idea, we consider Newton iteration i+l of the 
full-scale time evolution problem corresponding to parameter fi G V s at time step t n . The following 
solutions are supposed to be available: (i) the exact full-scale solution to the tangent problem (ii) 
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the reduced solution without system approximation (iii) the solution to the reduced tangent problem 
using the system approximation. In the second and latter cases, the linearisation of the equilibrium is 
performed around the solution to the full-scale system in order to eliminate the effect of accumulation 
of errors with time. For subdomain e £ £ icd , the distance between (i) and (iii) can be decomposed 
using a triangle inequality as follows: 



< 
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(67) 



Recall that V| e ' is the full-scale solution corresponding to the interior degrees of freedom of subdomain 
e i Cj 6 '' 6af^ is the interior solution obtained using the Galerkin POD model (no system approxima- 
tion), and 5a[ e ' is the solution obtained when employing the gappy technique. 

The value of the second term of the right-hand side of inequality (|68[) can be controlled by the rank 
of the "static" POD basis . It is reduced to zero in the limit case, where its rank is equal to the 

cardinality of set F s ^ e \ We choose the smallest rank of Dj 6 - 1 under the constraint that the following 
criterion holds for any iteration of the Newton algorithm, at any time step and for any parameter value 
in V s : 
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This technique ensures that the error due to the approximation of the nonlinear vectors of internal 
forces by the gappy technique is comparable to the one generated by the POD approximation of the 
displacements, at least for all parameter values in V s . 

Notice that this method requires the solution of the evolution problems corresponding to parameters 
in V s a certain number of times, but at reduced costs using the gappy technique. 



5.2.3 Selection of the control points 

For each subdomain e G £ rcd , we now have to decide which subset of interior nodes will be chosen as 
control points. This defines boolean operator P.| e \ and, together with obtained in the previous 

subsection, the required gappy operator G^K In the context of the DEIM, the selection is performed 

in a greedy manner, for increasing rank of operator Dj e ) , where we recall that the columns of this 
operator are hierarchically ordered by SVD. More precisely, at iteration j > of the greedy algorithm, 
the degree of freedom for which the gappy interpolation error 



is maximum is defined as a "control degree of freedom". Operator Dj 6 -''-' is composed of the j first 
columns of D^, while D^: +1 is the j + 1 th column of Dj 6 - 1 . Interpolation coefficient ft is obtained 



by solving the following optimisation problem: 



f3 3 = argmin 



dM^-BSVi pWJ> (71) 



The rank of the j th greedy iterate P.^ e ' J is j-times the number of scalar unknowns per interior node 
of subdomain e. In our implementation of the method, the node carrying the new "control degree of 
freedom" is added as a new "control point" , and all its associated degrees of freedom are controlled, 
which means that the corresponding entries in pj e )> J+1 are set to one. For an arbitrary subdomain e, 

the application of this method provides a number of "control nodes" equal to the rank of D: 6 -' . We 
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refer to reference [49| for more details about this technique, and in particular for a discussion about 
its optimality (in a greedy sense). 

The reduced integration domain obtained by applying the methodology described previously is 
represented in figure [7] and is used in the following result section. 

6 Results 

6.1 Observed speed-up 

We now solve the parametric, time-dependent lattice problem described in section 13^41 using the par- 
titioned model order reduction approach, and report the speed-up in terms of run time. We propose 
four different lattice structures, using 121 (figure |HJ, 256, 441 and 961 (figure [9} nodes for each of the 
10 subdomains. The snapshot that us used to compute the local reduced spaces is the one chosen 
in section 14.31 Let us recall that the cross-validation procedure leads us to omit any reduction in 
subdomain 6, whose associated subproblem will be solved exactly. The remainder of the subproblems 
are projected in the appropriate reduced spaces identified in section 13.41 using the Petrov-Galerkin 
formulation (system approximation) developed in section [5j We present speed-up results for the sim- 
ulations corresponding to 8 = 40° and 9 = 27°. These time solutions are not in the snapshot, and we 
can reasonably extrapolate that the observed speed-ups are representative of what can be expected for 
an arbitrary value of the parameter. 




Figure 8: Solution of the last time step of the fully discrete time-dependent problem for a load angle 
of 45°. The lattice structure represented here is composed of 121 nodes per subdomain. The darkest 
bars correspond to a completely damaged state of the material, while the lightest bars are undamaged. 

The proposed methodology is implemented in the commercial package Matlab, in a pseudo parallel 
fashion: the required operations that are local per subdomains are performed sequentially using a 
single processor. In this setting, we choose to solve the non-symmetric condensed interface problems 
using a direct LU factorisation. The reason for this is that no reduction of this problem has been 
developed so far. The number of interface degrees of freedom remains unchanged after the projection 
of the subproblems in the local reduced spaces. We therefore chose the implementation of the method 
that favored the observed speed-up, keeping in mind that it is pseudo-parallel. We will come back to 
this point in the conclusion of this work. 

To show the performance of the reduced order model, we first compute the fine solution of the fully 
discrete problem, without reduction, that corresponds to the first of the two particular load angles 
mentioned previously. The convergence tolerance for the Newton algorithm used at each time step 
(norm of the residual divided by the norm of the vector of external forces) is set to 10~ 7 . This is 
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Figure 9: Solution of the last time step of the fully discrete time-dependent problem for a load angle 
of 45°, using 961 lattice nodes per subdomain. 



the reference solution U cx . The accuracy of any approximate solution U app to this problem will be 
quantified using the following normalised error function: 



£ ||u app (^)-u ox (t„, M )| 



j,app,0) (u app ) 2 = . (72) 

t n eT h 

Secondly, an approximate solution U lnex is obtained by a straightforward time-reduction technique: 
the Newton algorithms are solved to a loose tolerance, and the error i/ a PP>(^)(U lncx ) is reported as a 
function of run time in figure [TO] This result is untitled "Full Scale Inexact" (notice that our use of 
the term "inexact" is not to be confused with the Inexact Newton Method, whereby one loosens the 
convergence tolerance of an iterative linear solver associated with the successive predictors of a Newton 
algorithm (6ll |) . 

Finally, we compare the speed-up obtained when using this straightforward approach to the one 
obtained with the projection-based partitioned reduction approaches. The error between the reference 
solution U cx and the one obtained by the Galerkin projection-based partitioned model order reduc- 
tion (no system approximation), denoted by U spod ; is the output v^v-M (U spod ) of the previously 
defined error function. The corresponding result is labelled "Partitioned POD" in figure [TU] The error 
;,app,(M) Qj pspod ) f solution U pspod obtained with the partitioned reduction technique and the system 
approximation is reported next, under the label "Partitioned POD + System Approximation". All 
these curves are reproduced for the second test load angle in figure [TT] 

The errors described previously are plotted for different levels of convergence of the Newton algo- 
rithms, in both the approximate full-scale case and the reduced cases, which provides a fair comparison 
ground for the various domain decomposition algorithms. 

Observing the two figures of results, the following remarks can be made: 

• a significant speed-up is obtained when using the partitioned model order reduction approach 
in the Petrov-Galerkin setting. This observation is only valid for certain range of accuracy. 
Indeed, the projection-based approach is limited, in terms of reachable accuracy, by the snapshot 
approximation of the POD transform, and by its subsequent truncation at a low order. For 
instance, in the top- right result of figure ITOl the error obtained with the reduction method cannot 
decrease under 2 x 10~ 3 . This is of course to be expected, and the remedy to this problem, if 
actually considered a problem, is to increase the size of the local reduced spaces. On the contrary, 
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(b) Relative error for the different models using 256 nodes 
per subdomain 




-"-Partitioned POD + SA •> 
— B — Partitioned POD * e 
- O - Full Scale Inexact * 

10 1 1 1 1 1 1 5 1 

100 200 300 400 500 600 700 

runtime (seconds) 

(d) Relative error for the different models using 961 nodes 
per subdomain 



Figure 10: Relative error for the full-scale solver and reduced order models as a function of the runtime 
for a load angle 9 — 40°. The different points of the curves are generated by loosening the convergence 
of the Newton algorithms. 
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Figure 11: Relative error for the full-scale solver and reduced order models as a function of the runtime 
for a load angle 9 = 27° 
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the error versus CPU time corresponding to the full-scale system can reach machine precision 
when decreasing the convergence tolerance of the Newton solvers. 

• the Galerkin version of the partitioned POD approach produces insignificant speed-ups. This is a 
well-known fact. The number of degrees of freedom is reduced compared to the full-scale system, 
but the costly integrations of the reduced generalised forces over the spatial domain forbids any 
benefit in terms of computational gain over the full scale model. 

• the speed-up, observed in the valid region of reachable accuracy for the POD-based reduced 
order models, increases with the number of degrees of freedom of the full-scale problem. This 
can be easily explained. The cost of solving the full-scale problem increases with the number 
of fine-scale degrees of freedom. However, the dimension of the local reduced spaces does not 
depend on this fine-scale refinement, but on the statistical properties of the parametric problem. 
Typically, one would expect that the numerical cost associated with the reduction technique does 
not increase with the number of degrees of freedom of the full-scale solvers. In practice, this is 
not the case as some computational overhead penalises our implementation of the partitioned 
model order reduction approach, not the least of which is the fact that the condensed interface 
problem is not reduced. This will be discussed in the conclusion of the paper. 

Notice that in practice, the simulations using the reduced models with system approximation are 
only performed with the lowest tolerance threshold for the Newton algorithm. The intermediate run 
times have only been given for demonstration purposes. 

6.2 Remarks about the numerical efficiency of the system approximation 

We now present the previous speed-up results in a different form. The aim is to show the trend 
in computational gain as a function of the number of degrees of freedom of the full-scale problem, 
when using the proposed reduction approach, in a unique graph. In order to so, the speed-up results 
reported previously are reported in figure [12] as a function of the ratio between the number of elements 
of the lattice and the number of degrees of elements that are connected to control nodes of the system 
approximation. This ratio increases in a roughly linear manner with the number of degrees of freedom 
of the full-scale problem. The different point of the curve are the one obtained with the lattice models 
comprising respectively 64, 121, 256, 441 and 961 nodes per subdomain, with an appropriately low 
tolerance for the nonlinear solution algorithm. 

The increase in the speed-up as function of the number of degrees of freedom of the full-scale 
problem appears clearly in this form. But more importantly, the graph shows that the observed speed- 
up is directly related to the size of the reduced integration domain. As mentioned previously, this 
is a clear indication that the main factor that prevents us from obtaining further speed-up with the 
proposed method is the fact that the interface problem is not reduced, which requires to perform 
integrations over a large number of elements. This is a path to explore in order to bring the idea of 
reduced order modelling in a partitioned framework to its full capability in the context of fracture. 
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Figure 12: Evolution of the speed-up with the ratio of the number of elements in the structure over 
the number of elements comprising the reduced integration domain. 



7 Conclusion and perspectives 

In this paper, we have proposed a partitioned model order reduction strategy for parametrised problems 
of nonlinear fracture mechanics. The domain coupling has been performed using the tried and tested 
primal Schur-complement domain decomposition method. The local subproblems have been reduced 
by projection in low-dimensional subspaces obtained by the snapshot POD. We have shown that 
this approach permits to reduce, in a flexible manner, the computational cost associated with highly 
nonlinear problems. In particular: 

• the local reduced spaces are generated independently, and have independent dimensions, which 
allows us to focus the numerical effort where it is most needed. In fracture mechanics, subdomains 
that are close to highly damaged zones need a richer model to account for the effect of topological 
changes. The local POD transforms automatically generate local reduced spaces of relatively 
large dimensions in these zones. 

• the domain decomposition framework enables us to switch from reduced local solvers to full local 
solvers in a transparent manner. This is particularly useful for the subdomains that contain 
process zones, as a solution obtained by projection would be more expensive than a direct 
solution for a desirable accuracy. 

• the transitition between "offline" and "online" computations becomes flexible. The reduced 
models can be used in the zones where the local reduced spaces converge quickly when enriching 
the snapshot space, while still computing snapshots and refining the reduced models via a direct 
local solver in the remaining subdomains. 

We have shown that such a flexibility results in a significant speed-up in the case of parametric fracture 
mechanics problems. This speed-up naturally increases when the size of the highly damaged zone, in 
which the information is highly uncorrelated, is small compared to the scale of the structure. 
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This work is a step towards an optimal cost-reduction strategy for fracture. Further work needs to 
be done to increase the understanding, robustness and performance of the method. Two main research 
avenues are particularly interesting from our point of view. Firstly, the interface problem itself was not 
reduced, in our case to guarantee the interface kinematic compatibility. This results in a suboptimal 
reduced order model and, in the case of parallel computing, would generate expensive communications 
through the network. A reduction of the interface problem using the POD can be done but is neither 
elegant nor technically easy. Using a dual Schur-complement domain decomposition method would 
allow the kinematic approximation of the subproblems to include the interface as well. However, the 
difficulty would be rejected on the necessary reduction of the interface Lagrange multiplier space. This 
issue is our current direction of research. 

An other difficulty is the load balancing mismatch that would occur when using such a strategy in 
parallel. CPUs which support domains that are not reduced, or domains for which the corresponding 
subproblems need to be projected in a space of relatively high dimension, would require to perform 
mode operations. Hence, the domain partitioning itself should be performed jointly with the model 
reduction in order to distribute the load evenly. 

Finally, we outlined throughout the paper some points that need further investigations but which 
are not directly related to the topic of reduced model partitioning addressed in this paper. The optimal 
choice of the snapshot samples used to construct a posteriori reduced order models is currently a very 
active research area (see for instance the review [45j concerning the reduced basis method, or the new 
developments proposed in 5JJ in the case of the snapshot POD). For arbitrary type of nonlinearity, 
a clear answer to this problem is, to date, not available. We have used a technique based on cross- 
validation, which, admittedly, requires a decently fine snapshot space in order to provide a relevant 
error estimate. In addition, our technique does not help find particular zones of non-smoothness in the 
parameter space. It only provides a general trend for the projection error. Furthermore, an important 
point related to this issue is that the error criterion that has been used in this work are all based on 
global euclidean norms, without consideration for the physical phenomenon of interest. We believe 
that developing a "goal-oriented" domain-decomposition-based reduced order modelling would help 
alleviate a certain number of issues related to the certification of reduced order modelling for general 
nonlinearities. The same remark applies to system approximations. The community seems to converge 
towards gappy-type evaluations of the nonlinear terms of the balance equations. Schemes based on 
direct collocation have been less investigated than gappy reconstruction when applied to reduced basis 
approaches, but may lead to good performances. 
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